package net.wesc.ac.server;

/**
 * Debye-Huckel模型，由P. Debye和E. Huckel于1923年提出，模型把溶剂当作连续介质，只考虑离子间的
 * 静电相互作用，通过求解Possion-Boltzmann方程求得。
 * 模型方程参见WESC Activity Coefficient用户手册方程(4.1.1)
 * 模型中只有一个参数a，其物理意义是离子在溶液中的半径
 * A与温度的关系为A正比于T的-3/2次方
 * B与温度的关系为B正比于T的-1/2次方
 */
public class LMDH implements LMfunc
{  
  static final int a0 = 0; //radii
  static final double A0 = 0.509;
  static final double B0 = 3.291;
  static final double temp0 = 298.15;
  static int Zc;
  static int Za;
  static int Va;
  static int Vc;
  static int V;
  static double temp = temp0;
  static double A = A0;
  static double B = B0;
  
  public LMDH(int tZc, int tZa, double ttemp){
	  Zc = tZc;
	  Za = tZa;
	  Vc = Za;
	  Va = Zc;
	  V = Va+Vc;
	  temp = ttemp;
	  A = A0*Math.pow(temp0/temp, 1.5);
	  B = B0*Math.pow(temp0/temp, 0.5);
  }
  
  public double[] initial()
  {
    double[] a = new double[1];
    a[a0] = 0.3;
    return a;
  } //initial

  public double val(double[] x, double[] a){
	  double I = 0.5*(Va*x[0]*Math.pow(Za, 2)+Vc*x[0]*Math.pow(Zc, 2)); //I ionic strength
      return (-A*Math.abs(Za*Zc)*Math.sqrt(I))/(1+B*a[a0]*Math.sqrt(I));
  } //val
  
  public double grad(double[] x, double[] a, int a_k){
	  if (a_k == a0){
		  double I = 0.5*(Va*x[0]*Math.pow(Za, 2)+Vc*x[0]*Math.pow(Zc, 2)); //I ionic strength
		  return (A*I*B*Math.abs(Za*Zc))/(Math.pow(1+B*a[a0]*Math.sqrt(I),2));
	  }
	  else {
		  assert false;
		  return 0.;
	  }
  } //grad
} //DHTest
